{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Paddle训练PINNs指导二维NACA翼型的网格自适应细化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## **一、项目介绍**\n",
    "\n",
    "### 1. 相关背景\n",
    "计算流体力学（CFD）是21世纪流体力学领域的重要技术之一，使用数值方法在计算机中对流体力学的控制方程进行求解，从而可预测流场的流动。目前有多种成熟的CFD软件，比如FLUENT、Star-CD、OpenFoam、CFX等。为了求解流体力学的控制方程，一个必要的前提是有适合数值方法的离散网格。生成离散网格的软件也有多种，比如ICEM、GAMBIT、ANSYS、Star-CD等。离散网格的质量（Quality）是影响CFD求解的关键因素。但是在传统的CFD流程中，离散网格的生成仅有一些指导经验。为了生成一个高质量的离散网格，往往需要对CFD求解结果进行网格敏感性分析（Mesh Sensitiviy Study）。这是一个耗时且无法跳过的过程。\n",
    "\n",
    "近年来基于物理信息约束神经元网络（PINNs）的方法在流体力学领域较为流行。PINN的特色是在神经元网络的损失函数（Loss Function）里包含流体力学的控制方程，使得神经元网络可以代替CFD的求解过程。对比CFD软件，PINNs的求解速度较快，一部分是因为PINNs对CFD流程的简化，另一部分是因为PINNs的计算是在GPU上进行。但是，目前PINN求解精度还有待提高。\n",
    "\n",
    "### 2. 功能目标\n",
    "本课题的目的是基于PINNs进行CFD离散网格优化，即并不采用PINN代替高精度的CFD求解器，而是用PINNs的求解速度快以及其具备一定的求解精度的特点，加快寻找最优的离散网格，从而缩短整体CFD流程的速度。本质上，本课题是为了验证一个假说，即高质量的离散网格，也对应相对高精度的PINN。\n",
    "\n",
    "### 3. 意义\n",
    "本课题的意义是探索区别于传统方式生成网格的技术，验证PINNs进行CFD离散网格优化的可行性，为工业应用计算软件提供一个基于机器学习的高质量网格生成与细化的新路径。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 二、设计方案与实现流程\n",
    "\n",
    "### 1. 网格优化总体方案\n",
    "\n",
    "我们所提出的基于PINNs对网格进行优化的具体方案示意如图所示。Level表示优化的层级，Points/Nodes表示训练点集或者网格的节点，Mesh代表三角形或者四面体网格，CFD solutions为传统CFD得到的参考解。注意，这里的CFD solution是否调用是可自由选择的，故在图中使用虚线来表示，可以根据实际训练情况灵活地选择是否使用CFD数据来进行无监督学习（无CFD数据)或监督学习（有CFD数据)。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/0e1e6e2ea930489b9b2a2af2cf6c2be77530ff2a4323426cb0b157c898e93e56)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. 环境配置\n",
    "\n",
    "1. 由于该项目需要进行PINNs模型的训练，因此我们选择使用以paddle为后端的科学计算库DeepXDE。\n",
    "2. 此外，由于需要对一些网格文件（如.msh或.vtk）进行读取和写入，需要安装meshio。\n",
    "3. 从二维散点集通过Delauney或Constrained Delauney三角化生成非结构网格，我们选择使用Triangle库。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# 机器学习库(自动识别paddle作为后端)\n",
    "!pip install deepxde\n",
    "# 网格文件读写库\n",
    "!pip install meshio\n",
    "# 二维Delauney三角化生成库\n",
    "!pip install triangle\n",
    "# 更新一些必要的库\n",
    "!pip install --upgrade numpy\n",
    "!pip install matplotlib==3.5.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "import deepxde as dde\n",
    "import meshio\n",
    "import triangle as tr\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. 读取已有网格文件\n",
    "\n",
    "官方提供的网格文件naca001065_aoa0_box12_4.msh和naca633418_aoa0_box12_4.msh无法通过python库meshio读取，同时也无法直接导入Fluent或workbench中进行使用。但可以使用gmsh软件打开。因此，我们先将原网格文件通过gmsh转化为.bdf格式，导入workbench后设定好边界信息再导出为.msh格式的文件。此时的.msh网格文件可被meshio、gmsh、fluent、workbench等通用软件读取。接下来，我们通过`meshio`来实现从workbench导出的网格文件（.msh）的节点坐标读取，并将网格点进行分类：翼型点、内部点、远场点。最后保存为后续模型训练的残差点集。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "from work.read_mesh_tools import convert_and_remove_duplicates, find_internal_points, plot_points_xy\n",
    "# 读取.msh文件\n",
    "mesh = meshio.read('work/NACA0010-65/level1/0_Mesh_files/NACA001065_mesh(workbench).msh', file_format='ansys')\n",
    "# 获取cells单元的每个edge对应的节点索引\n",
    "cells = mesh.cells\n",
    "all_edge_data = cells[0].data\n",
    "farfield_edge_data = np.vstack((cells[1].data, cells[2].data, cells[3].data, cells[4].data))\n",
    "airfoil_edge_data = cells[5].data\n",
    "# 去掉索引中重复的元素，只保留节点的索引\n",
    "all_points_index = convert_and_remove_duplicates(all_edge_data)\n",
    "farfield_points_index = convert_and_remove_duplicates(farfield_edge_data)\n",
    "airfoil_points_index = convert_and_remove_duplicates(airfoil_edge_data)\n",
    "# 通过索引找到每个区域包含的节点坐标\n",
    "all_points1 = mesh.points[all_points_index]\n",
    "farfield_points1 = mesh.points[farfield_points_index]\n",
    "airfoil_points1 = mesh.points[airfoil_points_index]\n",
    "internal_points1 = find_internal_points(all_points1, farfield_points1, airfoil_points1)\n",
    "# Save the array as npy\n",
    "np.save('work/NACA0010-65/level1/0_Mesh_files/all_mesh_points.npy', all_points1)\n",
    "np.save('work/NACA0010-65/level1/0_Mesh_files/ordered_airfoil_points.npy', airfoil_points1)\n",
    "np.save('work/NACA0010-65/level1/0_Mesh_files/farfield_points.npy', farfield_points1)\n",
    "np.save('work/NACA0010-65/level1/0_Mesh_files/internal_points.npy', internal_points1)\n",
    "# 画出网格节点分布\n",
    "plot_points_xy(all_points1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**翼型Naca0010-65提取的网格节点共4511个，散点图如下：**\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/80b425f867394485b9fdf8bc9ac4584b1494e69fe93248d5a8c49c9673590960)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "from work.read_mesh_tools import convert_and_remove_duplicates, find_internal_points, plot_points_xy, reorder_points\n",
    "# 读取.msh文件\n",
    "mesh = meshio.read('work/NACA6334-18/level1/0_Mesh_files/NACA633418_aoa0_box12_4.msh', file_format='ansys')\n",
    "# 获取cells单元的每个edge对应的节点索引\n",
    "cells = mesh.cells\n",
    "all_edge_data = cells[0].data\n",
    "farfield_edge_data = np.vstack((cells[1].data, cells[2].data, cells[3].data, cells[4].data))\n",
    "airfoil_edge_data = cells[5].data\n",
    "# 去掉索引中重复的元素，只保留节点的索引\n",
    "all_points_index = convert_and_remove_duplicates(all_edge_data)\n",
    "farfield_points_index = convert_and_remove_duplicates(farfield_edge_data)\n",
    "airfoil_points_index = convert_and_remove_duplicates(airfoil_edge_data)\n",
    "# 通过索引找到每个区域包含的节点坐标\n",
    "all_points2 = mesh.points[all_points_index]\n",
    "farfield_points2 = mesh.points[farfield_points_index]\n",
    "original_airfoil_points2 = mesh.points[airfoil_points_index]\n",
    "internal_points2 = find_internal_points(all_points2, farfield_points2, original_airfoil_points2)\n",
    "# 对该翼型和四周的点进行排序\n",
    "airfoil_points2 = np.vstack((original_airfoil_points2[1:104,:], original_airfoil_points2[0:1,:], original_airfoil_points2[104:,:]))\n",
    "farfield_points2 = reorder_points(farfield_points2)\n",
    "# Save the array as npy\n",
    "np.save('work/NACA6334-18/level1/0_Mesh_files/all_mesh_points.npy', all_points2)\n",
    "np.save('work/NACA6334-18/level1/0_Mesh_files/ordered_airfoil_points.npy', airfoil_points2)\n",
    "np.save('work/NACA6334-18/level1/0_Mesh_files/farfield_points.npy', farfield_points2)\n",
    "np.save('work/NACA6334-18/level1/0_Mesh_files/internal_points.npy', internal_points2)\n",
    "# 画出网格节点分布\n",
    "plot_points_xy(all_points2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**翼型Naca6334-18提取的网格节点共4522个，散点图如下：**\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/fdde835b2f2c40a2bcfc6d0525d85563f0db6224136b4ae3b704d6507891aa7c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. 基于节点的三角形非结构网格生成\n",
    "\n",
    "对于二维网格点集，可使用Delaunay三角化将其转换为三角形非结构网格。对于包含翼型的网格，难点在于怎样使用Constrained Delaunay Triangulation算法，即考虑翼型边界约束的三角化（翼型内部不需要生成网格）。**后续所有的从坐标点生成三角形网格、或者对原网格进行加密的步骤均基于此算法。**\n",
    "\n",
    "下面我们使用python的非结构网格生成库Triangle来实现上一步提取的4511个网格节点的Delaunay三角化。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 无边界约束的Delaunay三角化：\n",
    "\n",
    "若直接使用triangle对所提取网格节点进行Delauney三角化生成网格，则也会在翼型内部生成不需要的三角形网格。\n",
    "代码和生成结果如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "from work.tri_plot_tools import plot_mesh\n",
    "all_points1 = np.load('work/NACA0010-65/level1/0_Mesh_files/all_mesh_points.npy')\n",
    "# Delaunay Triangulation\n",
    "A = dict(vertices = all_points1)\n",
    "B = tr.triangulate(A)\n",
    "tr.comparev(plt, A, B, figsize=(60, 45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "翼型内部产生冗余网格\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e3c54c3ed1404f80a9e7d7c504c14797dfafe55d1755475aaa3f3ec7d55eb705)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "from work.tri_plot_tools import plot_mesh\n",
    "all_points2 = np.load('work/NACA6334-18/level1/0_Mesh_files/all_mesh_points.npy')\n",
    "# Delaunay Triangulation\n",
    "A = dict(vertices = all_points2)\n",
    "B = tr.triangulate(A)\n",
    "tr.comparev(plt, A, B, figsize=(60, 45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "翼型内部产生冗余网格\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/9aadf91066464e4e903fa1d45af2e11133bb788a3c1a4fd39966e11d17463746)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 边界约束（Constrained） Delaunay三角化：\n",
    "\n",
    "为了避免翼型内部也被三角化而产生不需要的网格，我们需要**对翼型表面（内边界）和矩形四边（外边界）设置边界约束**。即设置一组翼型上点与点按次序首尾相连的索引，作为segments加入到函数中，同时设置翼型内部的某个点（0.5，0）作为holes，即禁止在翼型的内部生成三角网格。\n",
    "代码和生成结果如下所示："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "polygon1 = np.concatenate((airfoil_points1, farfield_points1, internal_points1))\n",
    "# 设置约束边界的索引，翼型边界和矩形边界\n",
    "inner_airfoil_index1 = np.hstack((np.arange(0, airfoil_points1.shape[0], 1).reshape(airfoil_points1.shape[0], 1), \n",
    "                                 np.arange(1, airfoil_points1.shape[0]+1, 1).reshape(airfoil_points1.shape[0], 1)))\n",
    "inner_airfoil_index1[airfoil_points1.shape[0]-1, 1]=0\n",
    "outer_farfield_index1 = np.hstack((np.arange(0, farfield_points1.shape[0], 1).reshape(farfield_points1.shape[0], 1), \n",
    "                                  np.arange(1, farfield_points1.shape[0] +1, 1).reshape(farfield_points1.shape[0], 1)))\n",
    "outer_farfield_index1[farfield_points1.shape[0]-1, 1]=0\n",
    "outer_farfield_index1 = outer_farfield_index1 + inner_airfoil_index1.shape[0]\n",
    "inner_outer_index1 = np.vstack((inner_airfoil_index1, outer_farfield_index1))\n",
    "# Constrained Delaunay Triangulation\n",
    "A = dict(vertices = polygon1, segments = inner_outer_index1, holes=[[0.5, 0]])\n",
    "B = tr.triangulate(A,'p')\n",
    "tr.comparev(plt, A, B, figsize=(60, 45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**可以看到，使用带约束的Delaunay三角化还原生成的网格，在约束翼型边界后内部不生成冗余的网格。该网格与被提取节点的原网格完全一致（4511 points，8657 cells）。**\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/c4c163a05a3b4704886a997bf9cfb0ad8f21f587704048538af6c0c3d46ab3e0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "polygon2 = np.concatenate((airfoil_points2, farfield_points2, internal_points2))\n",
    "# 设置约束边界的索引，翼型边界和矩形边界\n",
    "inner_airfoil_index = np.hstack((np.arange(0, airfoil_points2.shape[0], 1).reshape(airfoil_points2.shape[0], 1), \n",
    "                                 np.arange(1, airfoil_points2.shape[0]+1, 1).reshape(airfoil_points2.shape[0], 1)))\n",
    "inner_airfoil_index[airfoil_points2.shape[0]-1, 1]=0\n",
    "outer_farfield_index = np.hstack((np.arange(0, farfield_points2.shape[0], 1).reshape(farfield_points2.shape[0], 1), \n",
    "                                  np.arange(1, farfield_points2.shape[0] +1, 1).reshape(farfield_points2.shape[0], 1)))\n",
    "outer_farfield_index[farfield_points2.shape[0]-1, 1]=0\n",
    "outer_farfield_index = outer_farfield_index + inner_airfoil_index.shape[0]\n",
    "inner_outer_index2 = np.vstack((inner_airfoil_index, outer_farfield_index))\n",
    "# Constrained Delaunay Triangulation\n",
    "A = dict(vertices = polygon2, segments = inner_outer_index2, holes=[[0.5, 0]])\n",
    "B = tr.triangulate(A,'p')\n",
    "tr.comparev(plt, A, B, figsize=(60, 45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**可以看到，使用带约束的Delaunay三角化还原生成的网格，在约束翼型边界后内部不生成冗余的网格。该网格与被提取节点的原网格完全一致（4522 points，8676 cells）。**\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/9ec565a0a3de4b7a989f4f629aaf2e0287d2af1dcc674c7c83c54e15fc36ae21)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. 传统CFD软件求解翼型附近跨声速流动\n",
    "\n",
    "翼型左方有一个0.7马赫数的均匀来流。起初我们分别使用**OpenFOAM**和**ANSYS Fluent**对翼型NACA0010-65进行了模拟。但OpenFOAM求解可压缩的两个求解器rhoSimpleFoam和rhoCentralFoam在调节各种参数设置后均得不到较好的结果，\n",
    "\n",
    "#### rhoSimpleFoam无法收敛，极易发散\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/508e732b966640b79cde78ce2b985a7b637ddccf513f4156bf23b814f5f13b92)\n",
    "\n",
    "#### rhoCentralFoam可以收敛，但有振荡\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/a345a269ff1749d487e58bd7e421c4307bb2d38ff06d45ddbfa483b8d081d093)\n",
    "\n",
    "#### Fluent收敛且无震荡\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/7b78e1f120734bf1bbc53d7276a466fdb67fe25b31f5472a892400dbefad36fa)\n",
    "\n",
    "#### Fluent模拟第二个翼型NACA6334-18也收敛且无震荡\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e7707f91923a4bd0b26206f83ffd6bf23047c229e96d4e11b147d90c1fc21544)\n",
    "\n",
    "\n",
    "因此，我们使用ANSYS Fluent来求解基于赛方提供网格（Level1）以及后续不同优化阶段的网格（Level2 - n）的流场，所用求解器以及相关设置均固定。\n",
    "Fluent求解得到的数据同时也可根据需要提供给机器学习模型进行监督学习。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在每个level的子文件夹**1_Fluent_results**内运行以下代码可**读取NACA0010-65和Naca6334-18的level1的Fluent求解结果**，并保存成供PINNs训练的.npy数组："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "%cd work/NACA0010-65/level1/1_Fluent_results/\n",
    "!python read_fluent_results.py\n",
    "%cd /home/aistudio/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "%cd work/NACA6334-18/level1/1_Fluent_results/\n",
    "!python read_fluent_results.py\n",
    "%cd /home/aistudio/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. 基于Paddle训练PINNs求解流场\n",
    "\n",
    "接下来，我们使用百度飞桨开发的机器学习库PaddlePaddle为后端的DeepXDE来搭建PINNs神经网络模型，并且使用来自Fluent得到的流场结果（速度$u$和$v$，压力$p$，密度$\\rho$）进行监督学习。DeepXDE提供了建立神经网络结构，设置边界条件等功能，并可以使用Adam和L-BFGS等经典优化器来进行训练。\n",
    "\n",
    "对于$Ma=0.7$的翼型气动问题，我们选择**二维Euler方程**作为控制方程。无粘极限下描述可压缩流动的质量、动量和能量守恒的二维稳态欧拉方程可以写成以下保守形式：\n",
    "\n",
    "$$\n",
    "\\frac{\\partial }{{\\partial x}}\\left[ {\\begin{array}{*{20}{c}}{\\rho u}\\\\{\\rho {u^2} + p}\\\\{\\rho uv}\\\\{u{\\rm(E+p)}}\\end{array}} \\right] + \\frac{\\partial }{{\\partial y}}\\left[ {\\begin{array}{*{20}{c}}{\\rho v}\\\\{\\rho uv}\\\\{\\rho {v^2} + p}\\\\{v{\\rm(E+p)}}\\end{array}} \\right] = 0\n",
    "$$\n",
    "\n",
    "其中，$\\rho$ 为密度，$p$ 为压力，$u$和$v$ 是速度，$E$ 是总能量。为了封闭方程组，还需考虑描述压力和能量关系的状态方程:\n",
    "\n",
    "$$\n",
    "p = (\\gamma  - 1)\\left( {E - \\frac{1}{2}\\rho \\left( {{u^2} + {v^2}} \\right)} \\right)\n",
    "$$\n",
    "\n",
    "其中，$\\gamma=1.4$ 为绝热指数。另外，考虑到Fluent模拟采用的是空气的物性条件，我们必须对其进行无量纲化。特征物理量定义为来流的物理量，即温度$T_0$=283.24 K，声速 $a_0=337.09m/s$，压力 $p_0=73048 pa$，密度 $\\rho_0 = 0.9kg/m^3$。\n",
    "\n",
    "将以下无量纲物理量代入Euler方程：$u'=u/a_0$，$v'=v/a_0$， $p'=p/p_0$，$\\rho'=\\rho/\\rho_0$\n",
    "\n",
    "进而得到**无量纲的Euler方程** （为方便表示，则去掉 '）：\n",
    "\n",
    "$$\\frac{\\partial }{{\\partial x}}\\left[ {\\begin{array}{*{20}{c}}\n",
    "{\\rho u}\\\\\n",
    "{\\gamma \\rho {u^2} + p}\\\\\n",
    "{\\gamma \\rho uv}\\\\\n",
    "{u{\\rm E}}\n",
    "\\end{array}} \\right] + \\frac{\\partial }{{\\partial y}}\\left[ {\\begin{array}{*{20}{c}}\n",
    "{\\rho v}\\\\\n",
    "{\\gamma \\rho uv}\\\\\n",
    "{\\gamma \\rho {v^2} + p}\\\\\n",
    "{v{\\rm E}}\n",
    "\\end{array}} \\right] = 0\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "该案例使用PINNs以监督学习的方式来求解流场，此处我们直接利用Fluent得到的流场数据（速度、密度、压力）作为监督数据，搭建神经网络进行训练。设置PDEs与数据的权重为1：10，使用Adam和L-BFGS优化器共同训练模型，并采用学习率逐步递减的分段训练策略。PINNs的训练可以运行各个网格优化level文件中的该程序来完成："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 训练模型（注意：时间较长，可跳过该块，直接加载训练好的模型）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "execution": {
     "iopub.execute_input": "2023-11-06T05:18:18.765678Z",
     "iopub.status.busy": "2023-11-06T05:18:18.764689Z",
     "iopub.status.idle": "2023-11-06T06:18:18.593081Z",
     "shell.execute_reply": "2023-11-06T06:18:18.592210Z",
     "shell.execute_reply.started": "2023-11-06T05:18:18.765621Z"
    },
    "jupyter": {
     "outputs_hidden": true
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/aistudio/work/NACA0010-65/level1/2_PINN_training\r\n",
      "Using backend: paddle\r\n",
      "Other available backends: tensorflow.compat.v1, tensorflow, pytorch, jax.\r\n",
      "paddle supports more examples now and is recommended.\r\n",
      " \r\n",
      "Set the default float type to float32\r\n",
      "Figure(1920x1440)\r\n",
      "Warning: CSGDifference.uniform_points not implemented. Use random_points instead.\r\n",
      "[4511, 4511, 4511, 4511]\r\n",
      "4511\r\n",
      "W1106 13:18:27.116514   731 gpu_resources.cc:85] Please NOTE: device: 0, GPU Compute Capability: 7.0, Driver API Version: 12.0, Runtime API Version: 11.2\r\n",
      "W1106 13:18:27.120384   731 gpu_resources.cc:115] device: 0, cuDNN Version: 8.2.\r\n",
      "Compiling model...\r\n",
      "'compile' took 0.000215 s\r\n",
      "\r\n",
      "Warning: epochs is deprecated and will be removed in a future version. Use iterations instead.\r\n",
      "Training model...\r\n",
      "\r\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\r\n",
      "0         [5.88e-04, 8.29e-03, 9.06e-03, 2.23e-03, 5.71e+00, 1.13e+00, 8.81e+00, 1.11e+01]    [6.34e-04, 2.76e-03, 2.00e-03, 1.73e-03, 5.71e+00, 1.13e+00, 8.81e+00, 1.11e+01]    []  \r\n",
      "1000      [7.22e-03, 6.15e-03, 3.39e-02, 1.27e-02, 1.14e-01, 4.08e-02, 1.78e-01, 6.13e-02]    [1.15e-03, 1.87e-03, 3.23e-03, 4.72e-03, 1.14e-01, 4.08e-02, 1.78e-01, 6.13e-02]    []  \r\n",
      "2000      [6.62e-03, 8.93e-03, 2.34e-02, 7.57e-03, 3.76e-02, 2.56e-02, 4.72e-02, 1.79e-02]    [1.13e-03, 1.15e-03, 3.38e-03, 2.60e-03, 3.76e-02, 2.56e-02, 4.72e-02, 1.79e-02]    []  \r\n",
      "3000      [6.89e-03, 6.98e-03, 1.39e-02, 3.97e-03, 3.19e-02, 1.82e-02, 3.10e-02, 1.46e-02]    [5.75e-04, 5.62e-04, 1.14e-03, 9.92e-04, 3.19e-02, 1.82e-02, 3.10e-02, 1.46e-02]    []  \r\n",
      "4000      [8.57e-03, 1.69e-02, 1.17e-02, 3.34e-02, 2.83e-02, 1.39e-02, 2.32e-02, 1.27e-02]    [8.35e-04, 1.34e-03, 8.55e-04, 2.65e-03, 2.83e-02, 1.39e-02, 2.32e-02, 1.27e-02]    []  \r\n",
      "5000      [5.16e-03, 7.56e-03, 7.51e-03, 1.12e-02, 1.94e-02, 1.01e-02, 1.58e-02, 8.76e-03]    [3.54e-04, 6.61e-04, 6.82e-04, 1.29e-03, 1.94e-02, 1.01e-02, 1.58e-02, 8.76e-03]    []  \r\n",
      "6000      [3.41e-03, 4.08e-03, 4.61e-03, 3.17e-03, 1.41e-02, 7.57e-03, 1.18e-02, 6.61e-03]    [2.64e-04, 6.62e-04, 6.29e-04, 7.97e-04, 1.41e-02, 7.57e-03, 1.18e-02, 6.61e-03]    []  \r\n",
      "7000      [3.01e-03, 3.98e-03, 3.44e-03, 4.68e-03, 1.13e-02, 6.13e-03, 9.36e-03, 5.23e-03]    [2.18e-04, 8.90e-04, 5.77e-04, 6.28e-04, 1.13e-02, 6.13e-03, 9.36e-03, 5.23e-03]    []  \r\n",
      "8000      [3.21e-03, 4.29e-03, 3.00e-03, 6.68e-03, 8.91e-03, 5.19e-03, 7.22e-03, 4.00e-03]    [2.16e-04, 1.07e-03, 5.48e-04, 6.11e-04, 8.91e-03, 5.19e-03, 7.22e-03, 4.00e-03]    []  \r\n",
      "9000      [2.54e-03, 3.78e-03, 1.93e-03, 2.80e-03, 7.22e-03, 4.62e-03, 6.00e-03, 3.28e-03]    [2.18e-04, 8.08e-04, 5.38e-04, 5.33e-04, 7.22e-03, 4.62e-03, 6.00e-03, 3.28e-03]    []  \r\n",
      "10000     [3.06e-03, 3.86e-03, 1.83e-03, 4.63e-03, 6.24e-03, 4.12e-03, 4.93e-03, 2.72e-03]    [2.50e-04, 6.70e-04, 4.88e-04, 4.93e-04, 6.24e-03, 4.12e-03, 4.93e-03, 2.72e-03]    []  \r\n",
      "11000     [2.88e-03, 3.31e-03, 1.43e-03, 1.89e-03, 5.54e-03, 3.88e-03, 4.29e-03, 2.39e-03]    [2.90e-04, 4.25e-04, 3.80e-04, 4.19e-04, 5.54e-03, 3.88e-03, 4.29e-03, 2.39e-03]    []  \r\n",
      "12000     [2.83e-03, 3.25e-03, 1.18e-03, 1.59e-03, 5.12e-03, 3.64e-03, 3.96e-03, 2.20e-03]    [2.85e-04, 3.77e-04, 3.43e-04, 3.58e-04, 5.12e-03, 3.64e-03, 3.96e-03, 2.20e-03]    []  \r\n",
      "13000     [3.22e-03, 5.44e-03, 1.41e-03, 8.97e-03, 4.94e-03, 3.59e-03, 3.99e-03, 2.20e-03]    [2.71e-04, 4.47e-04, 3.02e-04, 6.01e-04, 4.94e-03, 3.59e-03, 3.99e-03, 2.20e-03]    []  \r\n",
      "14000     [3.14e-03, 3.49e-03, 1.05e-03, 2.84e-03, 4.75e-03, 3.37e-03, 3.52e-03, 1.96e-03]    [3.21e-04, 2.93e-04, 2.47e-04, 4.48e-04, 4.75e-03, 3.37e-03, 3.52e-03, 1.96e-03]    []  \r\n",
      "15000     [2.86e-03, 3.25e-03, 7.86e-04, 1.49e-03, 4.45e-03, 3.07e-03, 3.41e-03, 1.87e-03]    [2.66e-04, 3.31e-04, 2.35e-04, 2.33e-04, 4.45e-03, 3.07e-03, 3.41e-03, 1.87e-03]    []  \r\n",
      "16000     [3.18e-03, 3.75e-03, 8.35e-04, 4.88e-03, 4.37e-03, 2.87e-03, 3.49e-03, 1.84e-03]    [2.72e-04, 4.03e-04, 2.59e-04, 3.50e-04, 4.37e-03, 2.87e-03, 3.49e-03, 1.84e-03]    []  \r\n",
      "17000     [2.77e-03, 3.32e-03, 6.48e-04, 9.64e-04, 4.19e-03, 2.91e-03, 3.17e-03, 1.74e-03]    [2.62e-04, 3.43e-04, 2.16e-04, 2.96e-04, 4.19e-03, 2.91e-03, 3.17e-03, 1.74e-03]    []  \r\n",
      "18000     [3.36e-03, 4.57e-03, 1.12e-03, 6.55e-03, 4.27e-03, 2.76e-03, 3.11e-03, 1.71e-03]    [2.94e-04, 3.31e-04, 2.28e-04, 3.76e-04, 4.27e-03, 2.76e-03, 3.11e-03, 1.71e-03]    []  \r\n",
      "19000     [2.73e-03, 3.22e-03, 5.77e-04, 1.28e-03, 3.97e-03, 2.59e-03, 2.99e-03, 1.63e-03]    [2.44e-04, 3.20e-04, 2.10e-04, 1.63e-04, 3.97e-03, 2.59e-03, 2.99e-03, 1.63e-03]    []  \r\n",
      "20000     [2.90e-03, 3.59e-03, 7.25e-04, 3.72e-03, 3.88e-03, 2.43e-03, 2.99e-03, 1.60e-03]    [2.45e-04, 3.57e-04, 2.30e-04, 2.45e-04, 3.88e-03, 2.43e-03, 2.99e-03, 1.60e-03]    []  \r\n",
      "\r\n",
      "Best model at step 19000:\r\n",
      "  train loss: 1.90e-02\r\n",
      "  test loss: 1.21e-02\r\n",
      "  test metric: []\r\n",
      "\r\n",
      "'train' took 447.093908 s\r\n",
      "\r\n",
      "Saving loss history to ./Train_history/loss.dat ...\r\n",
      "Saving training data to ./Train_history/train.dat ...\r\n",
      "Saving test data to ./Train_history/test.dat ...\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Compiling model...\r\n",
      "'compile' took 0.000264 s\r\n",
      "\r\n",
      "Warning: epochs is deprecated and will be removed in a future version. Use iterations instead.\r\n",
      "Training model...\r\n",
      "\r\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\r\n",
      "20000     [2.90e-03, 3.59e-03, 7.25e-04, 3.72e-03, 3.88e-03, 2.43e-03, 2.99e-03, 1.60e-03]    [2.45e-04, 3.57e-04, 2.30e-04, 2.45e-04, 3.88e-03, 2.43e-03, 2.99e-03, 1.60e-03]    []  \r\n",
      "21000     [2.69e-03, 3.10e-03, 5.22e-04, 8.86e-04, 3.81e-03, 2.44e-03, 2.81e-03, 1.53e-03]    [2.53e-04, 3.08e-04, 2.09e-04, 1.65e-04, 3.81e-03, 2.44e-03, 2.81e-03, 1.53e-03]    []  \r\n",
      "22000     [2.64e-03, 3.13e-03, 5.17e-04, 8.45e-04, 3.73e-03, 2.33e-03, 2.71e-03, 1.48e-03]    [2.50e-04, 3.15e-04, 2.05e-04, 1.70e-04, 3.73e-03, 2.33e-03, 2.71e-03, 1.48e-03]    []  \r\n",
      "23000     [2.71e-03, 2.93e-03, 5.08e-04, 1.11e-03, 3.67e-03, 2.16e-03, 2.62e-03, 1.43e-03]    [2.67e-04, 2.84e-04, 1.91e-04, 1.84e-04, 3.67e-03, 2.16e-03, 2.62e-03, 1.43e-03]    []  \r\n",
      "24000     [2.74e-03, 3.48e-03, 5.81e-04, 2.91e-03, 3.58e-03, 2.18e-03, 2.56e-03, 1.40e-03]    [2.48e-04, 3.26e-04, 1.73e-04, 2.34e-04, 3.58e-03, 2.18e-03, 2.56e-03, 1.40e-03]    []  \r\n",
      "25000     [2.57e-03, 3.03e-03, 5.37e-04, 1.13e-03, 3.52e-03, 1.98e-03, 2.45e-03, 1.35e-03]    [2.47e-04, 3.08e-04, 1.74e-04, 1.54e-04, 3.52e-03, 1.98e-03, 2.45e-03, 1.35e-03]    []  \r\n",
      "26000     [2.56e-03, 2.86e-03, 4.70e-04, 9.25e-04, 3.43e-03, 1.87e-03, 2.37e-03, 1.31e-03]    [2.60e-04, 2.87e-04, 1.55e-04, 1.68e-04, 3.43e-03, 1.87e-03, 2.37e-03, 1.31e-03]    []  \r\n",
      "27000     [2.44e-03, 2.94e-03, 4.55e-04, 8.25e-04, 3.35e-03, 1.80e-03, 2.31e-03, 1.28e-03]    [2.44e-04, 3.12e-04, 1.46e-04, 1.53e-04, 3.35e-03, 1.80e-03, 2.31e-03, 1.28e-03]    []  \r\n",
      "28000     [2.45e-03, 2.82e-03, 6.37e-04, 1.24e-03, 3.29e-03, 1.72e-03, 2.25e-03, 1.25e-03]    [2.54e-04, 2.88e-04, 1.38e-04, 1.57e-04, 3.29e-03, 1.72e-03, 2.25e-03, 1.25e-03]    []  \r\n",
      "29000     [2.40e-03, 2.71e-03, 4.32e-04, 1.02e-03, 3.23e-03, 1.59e-03, 2.19e-03, 1.22e-03]    [2.53e-04, 3.00e-04, 1.21e-04, 1.50e-04, 3.23e-03, 1.59e-03, 2.19e-03, 1.22e-03]    []  \r\n",
      "30000     [2.47e-03, 2.82e-03, 4.69e-04, 1.57e-03, 3.17e-03, 1.64e-03, 2.12e-03, 1.19e-03]    [2.68e-04, 2.84e-04, 1.18e-04, 2.66e-04, 3.17e-03, 1.64e-03, 2.12e-03, 1.19e-03]    []  \r\n",
      "31000     [2.23e-03, 2.86e-03, 4.56e-04, 1.30e-03, 3.13e-03, 1.46e-03, 2.11e-03, 1.18e-03]    [2.22e-04, 3.44e-04, 1.14e-04, 1.09e-04, 3.13e-03, 1.46e-03, 2.11e-03, 1.18e-03]    []  \r\n",
      "32000     [2.19e-03, 2.84e-03, 4.52e-04, 1.25e-03, 3.09e-03, 1.46e-03, 2.07e-03, 1.16e-03]    [2.22e-04, 3.44e-04, 1.12e-04, 1.43e-04, 3.09e-03, 1.46e-03, 2.07e-03, 1.16e-03]    []  \r\n",
      "33000     [2.22e-03, 2.51e-03, 4.46e-04, 9.32e-04, 3.02e-03, 1.37e-03, 1.99e-03, 1.13e-03]    [2.40e-04, 3.04e-04, 1.12e-04, 1.31e-04, 3.02e-03, 1.37e-03, 1.99e-03, 1.13e-03]    []  \r\n",
      "34000     [2.18e-03, 2.42e-03, 4.08e-04, 7.90e-04, 2.97e-03, 1.32e-03, 1.96e-03, 1.11e-03]    [2.39e-04, 3.01e-04, 1.13e-04, 1.40e-04, 2.97e-03, 1.32e-03, 1.96e-03, 1.11e-03]    []  \r\n",
      "35000     [2.28e-03, 2.44e-03, 5.34e-04, 1.85e-03, 2.94e-03, 1.28e-03, 1.94e-03, 1.10e-03]    [2.38e-04, 2.97e-04, 1.18e-04, 1.47e-04, 2.94e-03, 1.28e-03, 1.94e-03, 1.10e-03]    []  \r\n",
      "36000     [2.10e-03, 2.32e-03, 3.94e-04, 6.79e-04, 2.90e-03, 1.26e-03, 1.91e-03, 1.09e-03]    [2.30e-04, 2.92e-04, 1.18e-04, 1.38e-04, 2.90e-03, 1.26e-03, 1.91e-03, 1.09e-03]    []  \r\n",
      "37000     [2.43e-03, 2.21e-03, 4.32e-04, 2.48e-03, 2.94e-03, 1.25e-03, 1.89e-03, 1.09e-03]    [2.63e-04, 2.35e-04, 1.16e-04, 2.16e-04, 2.94e-03, 1.25e-03, 1.89e-03, 1.09e-03]    []  \r\n",
      "38000     [2.04e-03, 2.28e-03, 4.03e-04, 8.55e-04, 2.87e-03, 1.22e-03, 1.88e-03, 1.07e-03]    [2.16e-04, 2.84e-04, 1.20e-04, 1.17e-04, 2.87e-03, 1.22e-03, 1.88e-03, 1.07e-03]    []  \r\n",
      "39000     [2.09e-03, 2.94e-03, 7.03e-04, 3.00e-03, 2.82e-03, 1.19e-03, 1.93e-03, 1.07e-03]    [1.95e-04, 3.54e-04, 1.13e-04, 1.60e-04, 2.82e-03, 1.19e-03, 1.93e-03, 1.07e-03]    []  \r\n",
      "40000     [2.09e-03, 2.88e-03, 1.13e-03, 1.92e-03, 2.84e-03, 1.21e-03, 1.97e-03, 1.10e-03]    [2.15e-04, 3.34e-04, 1.52e-04, 2.22e-04, 2.84e-03, 1.21e-03, 1.97e-03, 1.10e-03]    []  \r\n",
      "\r\n",
      "Best model at step 38000:\r\n",
      "  train loss: 1.26e-02\r\n",
      "  test loss: 7.77e-03\r\n",
      "  test metric: []\r\n",
      "\r\n",
      "'train' took 488.360785 s\r\n",
      "\r\n",
      "Saving loss history to ./Train_history/loss.dat ...\r\n",
      "Saving training data to ./Train_history/train.dat ...\r\n",
      "Saving test data to ./Train_history/test.dat ...\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Compiling model...\r\n",
      "'compile' took 0.000355 s\r\n",
      "\r\n",
      "Warning: epochs is deprecated and will be removed in a future version. Use iterations instead.\r\n",
      "Training model...\r\n",
      "\r\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\r\n",
      "40000     [2.09e-03, 2.88e-03, 1.13e-03, 1.92e-03, 2.84e-03, 1.21e-03, 1.97e-03, 1.10e-03]    [2.15e-04, 3.34e-04, 1.52e-04, 2.22e-04, 2.84e-03, 1.21e-03, 1.97e-03, 1.10e-03]    []  \r\n",
      "41000     [2.00e-03, 2.14e-03, 3.62e-04, 5.67e-04, 2.77e-03, 1.15e-03, 1.82e-03, 1.04e-03]    [2.11e-04, 2.81e-04, 1.09e-04, 1.23e-04, 2.77e-03, 1.15e-03, 1.82e-03, 1.04e-03]    []  \r\n",
      "42000     [1.98e-03, 2.11e-03, 3.60e-04, 5.59e-04, 2.75e-03, 1.13e-03, 1.80e-03, 1.03e-03]    [2.07e-04, 2.81e-04, 1.08e-04, 1.22e-04, 2.75e-03, 1.13e-03, 1.80e-03, 1.03e-03]    []  \r\n",
      "43000     [1.96e-03, 2.08e-03, 3.57e-04, 5.60e-04, 2.73e-03, 1.11e-03, 1.79e-03, 1.02e-03]    [2.03e-04, 2.80e-04, 1.06e-04, 1.22e-04, 2.73e-03, 1.11e-03, 1.79e-03, 1.02e-03]    []  \r\n",
      "44000     [1.94e-03, 2.04e-03, 3.53e-04, 5.44e-04, 2.71e-03, 1.10e-03, 1.77e-03, 1.01e-03]    [1.99e-04, 2.75e-04, 1.05e-04, 1.18e-04, 2.71e-03, 1.10e-03, 1.77e-03, 1.01e-03]    []  \r\n",
      "45000     [1.91e-03, 2.03e-03, 3.50e-04, 5.33e-04, 2.68e-03, 1.09e-03, 1.76e-03, 1.00e-03]    [1.93e-04, 2.74e-04, 1.04e-04, 1.13e-04, 2.68e-03, 1.09e-03, 1.76e-03, 1.00e-03]    []  \r\n",
      "46000     [1.88e-03, 2.01e-03, 3.52e-04, 5.44e-04, 2.66e-03, 1.09e-03, 1.75e-03, 9.93e-04]    [1.87e-04, 2.68e-04, 1.05e-04, 1.08e-04, 2.66e-03, 1.09e-03, 1.75e-03, 9.93e-04]    []  \r\n",
      "47000     [1.86e-03, 1.95e-03, 3.45e-04, 5.24e-04, 2.63e-03, 1.07e-03, 1.73e-03, 9.84e-04]    [1.85e-04, 2.63e-04, 1.04e-04, 1.07e-04, 2.63e-03, 1.07e-03, 1.73e-03, 9.84e-04]    []  \r\n",
      "48000     [1.87e-03, 1.88e-03, 3.42e-04, 5.84e-04, 2.60e-03, 1.05e-03, 1.72e-03, 9.76e-04]    [1.87e-04, 2.52e-04, 1.02e-04, 1.11e-04, 2.60e-03, 1.05e-03, 1.72e-03, 9.76e-04]    []  \r\n",
      "49000     [1.80e-03, 1.89e-03, 3.41e-04, 5.09e-04, 2.58e-03, 1.05e-03, 1.71e-03, 9.66e-04]    [1.77e-04, 2.55e-04, 1.03e-04, 9.97e-05, 2.58e-03, 1.05e-03, 1.71e-03, 9.66e-04]    []  \r\n",
      "50000     [1.70e-03, 2.09e-03, 3.54e-04, 7.22e-04, 2.55e-03, 1.04e-03, 1.72e-03, 9.61e-04]    [1.57e-04, 2.78e-04, 1.05e-04, 9.73e-05, 2.55e-03, 1.04e-03, 1.72e-03, 9.61e-04]    []  \r\n",
      "51000     [1.75e-03, 1.82e-03, 3.54e-04, 5.85e-04, 2.51e-03, 1.05e-03, 1.69e-03, 9.46e-04]    [1.67e-04, 2.42e-04, 1.04e-04, 9.58e-05, 2.51e-03, 1.05e-03, 1.69e-03, 9.46e-04]    []  \r\n",
      "52000     [1.70e-03, 1.75e-03, 3.57e-04, 5.56e-04, 2.47e-03, 1.04e-03, 1.67e-03, 9.35e-04]    [1.59e-04, 2.34e-04, 1.07e-04, 9.21e-05, 2.47e-03, 1.04e-03, 1.67e-03, 9.35e-04]    []  \r\n",
      "53000     [1.63e-03, 1.71e-03, 3.72e-04, 4.85e-04, 2.43e-03, 1.03e-03, 1.66e-03, 9.21e-04]    [1.47e-04, 2.34e-04, 1.10e-04, 9.67e-05, 2.43e-03, 1.03e-03, 1.66e-03, 9.21e-04]    []  \r\n",
      "54000     [1.55e-03, 1.63e-03, 3.93e-04, 4.77e-04, 2.38e-03, 1.02e-03, 1.63e-03, 9.04e-04]    [1.34e-04, 2.28e-04, 1.15e-04, 1.02e-04, 2.38e-03, 1.02e-03, 1.63e-03, 9.04e-04]    []  \r\n",
      "55000     [1.44e-03, 1.55e-03, 4.18e-04, 4.62e-04, 2.33e-03, 1.02e-03, 1.60e-03, 8.86e-04]    [1.22e-04, 2.25e-04, 1.20e-04, 1.16e-04, 2.33e-03, 1.02e-03, 1.60e-03, 8.86e-04]    []  \r\n",
      "56000     [1.35e-03, 1.45e-03, 4.41e-04, 4.74e-04, 2.28e-03, 1.02e-03, 1.57e-03, 8.67e-04]    [1.15e-04, 2.21e-04, 1.21e-04, 1.20e-04, 2.28e-03, 1.02e-03, 1.57e-03, 8.67e-04]    []  \r\n",
      "57000     [1.24e-03, 1.36e-03, 4.57e-04, 4.57e-04, 2.23e-03, 1.01e-03, 1.54e-03, 8.49e-04]    [1.11e-04, 2.18e-04, 1.18e-04, 1.40e-04, 2.23e-03, 1.01e-03, 1.54e-03, 8.49e-04]    []  \r\n",
      "58000     [1.15e-03, 1.27e-03, 4.61e-04, 4.30e-04, 2.19e-03, 1.00e-03, 1.51e-03, 8.35e-04]    [1.07e-04, 2.09e-04, 1.15e-04, 1.49e-04, 2.19e-03, 1.00e-03, 1.51e-03, 8.35e-04]    []  \r\n",
      "59000     [1.08e-03, 1.23e-03, 4.82e-04, 6.14e-04, 2.15e-03, 1.02e-03, 1.48e-03, 8.20e-04]    [9.94e-05, 1.99e-04, 1.14e-04, 1.28e-04, 2.15e-03, 1.02e-03, 1.48e-03, 8.20e-04]    []  \r\n",
      "60000     [1.02e-03, 1.15e-03, 4.76e-04, 5.68e-04, 2.13e-03, 9.90e-04, 1.47e-03, 8.09e-04]    [9.82e-05, 1.97e-04, 1.18e-04, 1.50e-04, 2.13e-03, 9.90e-04, 1.47e-03, 8.09e-04]    []  \r\n",
      "\r\n",
      "Best model at step 60000:\r\n",
      "  train loss: 8.61e-03\r\n",
      "  test loss: 5.95e-03\r\n",
      "  test metric: []\r\n",
      "\r\n",
      "'train' took 479.395846 s\r\n",
      "\r\n",
      "Saving loss history to ./Train_history/loss.dat ...\r\n",
      "Saving training data to ./Train_history/train.dat ...\r\n",
      "Saving test data to ./Train_history/test.dat ...\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Compiling model...\r\n",
      "'compile' took 0.000333 s\r\n",
      "\r\n",
      "Training model...\r\n",
      "\r\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\r\n",
      "60000     [1.02e-03, 1.15e-03, 4.76e-04, 5.68e-04, 2.13e-03, 9.90e-04, 1.47e-03, 8.09e-04]    [9.82e-05, 1.97e-04, 1.18e-04, 1.50e-04, 2.13e-03, 9.90e-04, 1.47e-03, 8.09e-04]    []  \r\n",
      "61000     [5.57e-04, 7.57e-04, 3.10e-04, 2.67e-04, 1.94e-03, 9.44e-04, 1.26e-03, 7.26e-04]    [6.02e-05, 1.26e-04, 2.22e-04, 2.58e-04, 1.94e-03, 9.44e-04, 1.26e-03, 7.26e-04]    []  \r\n",
      "Epoch 61000: train loss improved from inf to 6.76e-03, saving model to Model/model-61000.pdparams ...\r\n",
      "\r\n",
      "62000     [4.33e-04, 6.06e-04, 2.07e-04, 2.06e-04, 1.83e-03, 8.96e-04, 1.21e-03, 7.00e-04]    [6.02e-05, 9.90e-05, 8.46e-05, 1.74e-04, 1.83e-03, 8.96e-04, 1.21e-03, 7.00e-04]    []  \r\n",
      "Epoch 62000: train loss improved from 6.76e-03 to 6.09e-03, saving model to Model/model-62000.pdparams ...\r\n",
      "\r\n",
      "63000     [3.88e-04, 5.39e-04, 1.66e-04, 1.87e-04, 1.73e-03, 8.67e-04, 1.17e-03, 6.68e-04]    [7.56e-05, 1.17e-04, 1.73e-04, 2.82e-04, 1.73e-03, 8.67e-04, 1.17e-03, 6.68e-04]    []  \r\n",
      "Epoch 63000: train loss improved from 6.09e-03 to 5.72e-03, saving model to Model/model-63000.pdparams ...\r\n",
      "\r\n",
      "64000     [3.37e-04, 4.74e-04, 1.39e-04, 1.69e-04, 1.67e-03, 8.77e-04, 1.13e-03, 6.49e-04]    [7.64e-05, 1.32e-04, 1.98e-04, 3.43e-04, 1.67e-03, 8.77e-04, 1.13e-03, 6.49e-04]    []  \r\n",
      "Epoch 64000: train loss improved from 5.72e-03 to 5.44e-03, saving model to Model/model-64000.pdparams ...\r\n",
      "\r\n",
      "65000     [3.16e-04, 4.46e-04, 1.35e-04, 1.63e-04, 1.61e-03, 8.48e-04, 1.08e-03, 6.27e-04]    [1.51e-04, 1.38e-04, 4.33e-04, 7.48e-04, 1.61e-03, 8.48e-04, 1.08e-03, 6.27e-04]    []  \r\n",
      "Epoch 65000: train loss improved from 5.44e-03 to 5.22e-03, saving model to Model/model-65000.pdparams ...\r\n",
      "\r\n",
      "66000     [3.03e-04, 4.28e-04, 1.15e-04, 1.44e-04, 1.55e-03, 8.24e-04, 1.04e-03, 6.05e-04]    [1.59e-04, 1.55e-04, 7.85e-04, 7.26e-04, 1.55e-03, 8.24e-04, 1.04e-03, 6.05e-04]    []  \r\n",
      "Epoch 66000: train loss improved from 5.22e-03 to 5.01e-03, saving model to Model/model-66000.pdparams ...\r\n",
      "\r\n",
      "67000     [2.86e-04, 4.15e-04, 1.08e-04, 1.39e-04, 1.51e-03, 8.15e-04, 1.00e-03, 5.84e-04]    [1.42e-04, 1.69e-04, 1.61e-03, 7.23e-04, 1.51e-03, 8.15e-04, 1.00e-03, 5.84e-04]    []  \r\n",
      "Epoch 67000: train loss improved from 5.01e-03 to 4.85e-03, saving model to Model/model-67000.pdparams ...\r\n",
      "\r\n",
      "68000     [2.71e-04, 3.90e-04, 1.07e-04, 1.40e-04, 1.47e-03, 7.99e-04, 9.74e-04, 5.69e-04]    [1.60e-04, 2.73e-04, 1.67e-03, 9.00e-04, 1.47e-03, 7.99e-04, 9.74e-04, 5.69e-04]    []  \r\n",
      "Epoch 68000: train loss improved from 4.85e-03 to 4.72e-03, saving model to Model/model-68000.pdparams ...\r\n",
      "\r\n",
      "69000     [2.84e-04, 3.98e-04, 1.14e-04, 1.42e-04, 1.42e-03, 7.67e-04, 9.27e-04, 5.43e-04]    [1.30e-04, 2.66e-04, 1.62e-03, 6.66e-04, 1.42e-03, 7.67e-04, 9.27e-04, 5.43e-04]    []  \r\n",
      "Epoch 69000: train loss improved from 4.72e-03 to 4.60e-03, saving model to Model/model-69000.pdparams ...\r\n",
      "\r\n",
      "70000     [2.72e-04, 3.94e-04, 1.10e-04, 1.47e-04, 1.40e-03, 7.48e-04, 8.90e-04, 5.23e-04]    [1.05e-04, 2.73e-04, 1.59e-03, 5.79e-04, 1.40e-03, 7.48e-04, 8.90e-04, 5.23e-04]    []  \r\n",
      "Epoch 70000: train loss improved from 4.60e-03 to 4.48e-03, saving model to Model/model-70000.pdparams ...\r\n",
      "\r\n",
      "71000     [2.72e-04, 3.95e-04, 1.03e-04, 1.44e-04, 1.38e-03, 7.19e-04, 8.54e-04, 5.04e-04]    [1.21e-04, 2.35e-04, 1.44e-03, 5.96e-04, 1.38e-03, 7.19e-04, 8.54e-04, 5.04e-04]    []  \r\n",
      "Epoch 71000: train loss improved from 4.48e-03 to 4.37e-03, saving model to Model/model-71000.pdparams ...\r\n",
      "\r\n",
      "72000     [2.77e-04, 3.97e-04, 1.13e-04, 1.42e-04, 1.35e-03, 6.90e-04, 8.12e-04, 4.85e-04]    [9.82e-05, 2.02e-04, 1.59e-03, 4.50e-04, 1.35e-03, 6.90e-04, 8.12e-04, 4.85e-04]    []  \r\n",
      "Epoch 72000: train loss improved from 4.37e-03 to 4.27e-03, saving model to Model/model-72000.pdparams ...\r\n",
      "\r\n",
      "73000     [2.74e-04, 3.99e-04, 1.16e-04, 1.45e-04, 1.32e-03, 6.59e-04, 7.70e-04, 4.67e-04]    [9.57e-05, 2.54e-04, 1.81e-03, 4.75e-04, 1.32e-03, 6.59e-04, 7.70e-04, 4.67e-04]    []  \r\n",
      "Epoch 73000: train loss improved from 4.27e-03 to 4.15e-03, saving model to Model/model-73000.pdparams ...\r\n",
      "\r\n",
      "74000     [2.71e-04, 3.96e-04, 1.20e-04, 1.51e-04, 1.30e-03, 6.33e-04, 7.39e-04, 4.49e-04]    [1.21e-04, 3.24e-04, 1.79e-03, 6.71e-04, 1.30e-03, 6.33e-04, 7.39e-04, 4.49e-04]    []  \r\n",
      "Epoch 74000: train loss improved from 4.15e-03 to 4.06e-03, saving model to Model/model-74000.pdparams ...\r\n",
      "\r\n",
      "75000     [2.68e-04, 3.93e-04, 1.20e-04, 1.52e-04, 1.27e-03, 6.20e-04, 7.11e-04, 4.35e-04]    [1.18e-04, 2.93e-04, 1.68e-03, 5.99e-04, 1.27e-03, 6.20e-04, 7.11e-04, 4.35e-04]    []  \r\n",
      "Epoch 75000: train loss improved from 4.06e-03 to 3.96e-03, saving model to Model/model-75000.pdparams ...\r\n",
      "\r\n",
      "76000     [2.68e-04, 3.97e-04, 1.17e-04, 1.59e-04, 1.24e-03, 5.91e-04, 6.77e-04, 4.19e-04]    [1.41e-04, 2.88e-04, 1.26e-03, 6.77e-04, 1.24e-03, 5.91e-04, 6.77e-04, 4.19e-04]    []  \r\n",
      "Epoch 76000: train loss improved from 3.96e-03 to 3.87e-03, saving model to Model/model-76000.pdparams ...\r\n",
      "\r\n",
      "77000     [2.63e-04, 3.93e-04, 1.16e-04, 1.51e-04, 1.23e-03, 5.74e-04, 6.56e-04, 4.06e-04]    [1.90e-04, 2.90e-04, 9.59e-04, 9.92e-04, 1.23e-03, 5.74e-04, 6.56e-04, 4.06e-04]    []  \r\n",
      "Epoch 77000: train loss improved from 3.87e-03 to 3.79e-03, saving model to Model/model-77000.pdparams ...\r\n",
      "\r\n",
      "78000     [2.59e-04, 3.87e-04, 1.14e-04, 1.48e-04, 1.21e-03, 5.63e-04, 6.37e-04, 3.97e-04]    [2.30e-04, 3.10e-04, 9.41e-04, 1.24e-03, 1.21e-03, 5.63e-04, 6.37e-04, 3.97e-04]    []  \r\n",
      "Epoch 78000: train loss improved from 3.79e-03 to 3.71e-03, saving model to Model/model-78000.pdparams ...\r\n",
      "\r\n",
      "79000     [2.55e-04, 3.89e-04, 1.02e-04, 1.48e-04, 1.19e-03, 5.62e-04, 6.21e-04, 3.86e-04]    [1.90e-04, 3.15e-04, 9.61e-04, 9.68e-04, 1.19e-03, 5.62e-04, 6.21e-04, 3.86e-04]    []  \r\n",
      "Epoch 79000: train loss improved from 3.71e-03 to 3.65e-03, saving model to Model/model-79000.pdparams ...\r\n",
      "\r\n",
      "80000     [2.60e-04, 3.74e-04, 1.02e-04, 1.34e-04, 1.18e-03, 5.51e-04, 6.10e-04, 3.80e-04]    [2.27e-04, 3.73e-04, 1.08e-03, 1.17e-03, 1.18e-03, 5.51e-04, 6.10e-04, 3.80e-04]    []  \r\n",
      "Epoch 80000: train loss improved from 3.65e-03 to 3.59e-03, saving model to Model/model-80000.pdparams ...\r\n",
      "\r\n",
      "81000     [2.60e-04, 3.75e-04, 1.00e-04, 1.35e-04, 1.16e-03, 5.40e-04, 5.94e-04, 3.72e-04]    [3.27e-04, 4.30e-04, 1.02e-03, 1.77e-03, 1.16e-03, 5.40e-04, 5.94e-04, 3.72e-04]    []  \r\n",
      "Epoch 81000: train loss improved from 3.59e-03 to 3.54e-03, saving model to Model/model-81000.pdparams ...\r\n",
      "\r\n",
      "82000     [2.47e-04, 3.68e-04, 9.52e-05, 1.38e-04, 1.16e-03, 5.39e-04, 5.87e-04, 3.67e-04]    [3.57e-04, 4.68e-04, 8.95e-04, 1.94e-03, 1.16e-03, 5.39e-04, 5.87e-04, 3.67e-04]    []  \r\n",
      "Epoch 82000: train loss improved from 3.54e-03 to 3.50e-03, saving model to Model/model-82000.pdparams ...\r\n",
      "\r\n",
      "83000     [2.51e-04, 3.65e-04, 9.49e-05, 1.42e-04, 1.14e-03, 5.26e-04, 5.75e-04, 3.59e-04]    [3.55e-04, 4.94e-04, 1.12e-03, 1.80e-03, 1.14e-03, 5.26e-04, 5.75e-04, 3.59e-04]    []  \r\n",
      "Epoch 83000: train loss improved from 3.50e-03 to 3.46e-03, saving model to Model/model-83000.pdparams ...\r\n",
      "\r\n",
      "84000     [2.41e-04, 3.62e-04, 9.78e-05, 1.44e-04, 1.13e-03, 5.21e-04, 5.65e-04, 3.55e-04]    [3.14e-04, 4.51e-04, 1.26e-03, 1.46e-03, 1.13e-03, 5.21e-04, 5.65e-04, 3.55e-04]    []  \r\n",
      "Epoch 84000: train loss improved from 3.46e-03 to 3.42e-03, saving model to Model/model-84000.pdparams ...\r\n",
      "\r\n",
      "85000     [2.37e-04, 3.51e-04, 1.00e-04, 1.46e-04, 1.13e-03, 5.15e-04, 5.55e-04, 3.49e-04]    [2.46e-04, 4.21e-04, 1.36e-03, 1.12e-03, 1.13e-03, 5.15e-04, 5.55e-04, 3.49e-04]    []  \r\n",
      "Epoch 85000: train loss improved from 3.42e-03 to 3.38e-03, saving model to Model/model-85000.pdparams ...\r\n",
      "\r\n",
      "86000     [2.37e-04, 3.53e-04, 9.89e-05, 1.47e-04, 1.11e-03, 5.06e-04, 5.43e-04, 3.43e-04]    [2.13e-04, 3.99e-04, 1.26e-03, 9.62e-04, 1.11e-03, 5.06e-04, 5.43e-04, 3.43e-04]    []  \r\n",
      "Epoch 86000: train loss improved from 3.38e-03 to 3.34e-03, saving model to Model/model-86000.pdparams ...\r\n",
      "\r\n",
      "87000     [2.35e-04, 3.47e-04, 9.79e-05, 1.46e-04, 1.10e-03, 4.97e-04, 5.35e-04, 3.38e-04]    [1.94e-04, 3.77e-04, 1.42e-03, 8.21e-04, 1.10e-03, 4.97e-04, 5.35e-04, 3.38e-04]    []  \r\n",
      "Epoch 87000: train loss improved from 3.34e-03 to 3.30e-03, saving model to Model/model-87000.pdparams ...\r\n",
      "\r\n",
      "88000     [2.27e-04, 3.41e-04, 9.70e-05, 1.41e-04, 1.08e-03, 4.97e-04, 5.28e-04, 3.36e-04]    [1.79e-04, 4.26e-04, 1.60e-03, 8.34e-04, 1.08e-03, 4.97e-04, 5.28e-04, 3.36e-04]    []  \r\n",
      "Epoch 88000: train loss improved from 3.30e-03 to 3.25e-03, saving model to Model/model-88000.pdparams ...\r\n",
      "\r\n",
      "89000     [2.20e-04, 3.31e-04, 1.01e-04, 1.40e-04, 1.08e-03, 4.93e-04, 5.21e-04, 3.34e-04]    [1.97e-04, 5.00e-04, 1.84e-03, 1.01e-03, 1.08e-03, 4.93e-04, 5.21e-04, 3.34e-04]    []  \r\n",
      "Epoch 89000: train loss improved from 3.25e-03 to 3.22e-03, saving model to Model/model-89000.pdparams ...\r\n",
      "\r\n",
      "90000     [2.13e-04, 3.19e-04, 1.03e-04, 1.40e-04, 1.06e-03, 4.92e-04, 5.14e-04, 3.32e-04]    [2.30e-04, 6.28e-04, 1.87e-03, 1.25e-03, 1.06e-03, 4.92e-04, 5.14e-04, 3.32e-04]    []  \r\n",
      "Epoch 90000: train loss improved from 3.22e-03 to 3.18e-03, saving model to Model/model-90000.pdparams ...\r\n",
      "\r\n",
      "91000     [2.07e-04, 3.15e-04, 1.02e-04, 1.35e-04, 1.06e-03, 4.80e-04, 5.08e-04, 3.30e-04]    [2.47e-04, 7.24e-04, 2.51e-03, 1.45e-03, 1.06e-03, 4.80e-04, 5.08e-04, 3.30e-04]    []  \r\n",
      "Epoch 91000: train loss improved from 3.18e-03 to 3.13e-03, saving model to Model/model-91000.pdparams ...\r\n",
      "\r\n",
      "92000     [2.05e-04, 3.07e-04, 1.05e-04, 1.36e-04, 1.05e-03, 4.68e-04, 4.99e-04, 3.24e-04]    [2.85e-04, 7.51e-04, 2.95e-03, 1.63e-03, 1.05e-03, 4.68e-04, 4.99e-04, 3.24e-04]    []  \r\n",
      "Epoch 92000: train loss improved from 3.13e-03 to 3.09e-03, saving model to Model/model-92000.pdparams ...\r\n",
      "\r\n",
      "93000     [2.04e-04, 3.03e-04, 1.00e-04, 1.35e-04, 1.04e-03, 4.58e-04, 4.92e-04, 3.20e-04]    [2.90e-04, 7.68e-04, 3.16e-03, 1.54e-03, 1.04e-03, 4.58e-04, 4.92e-04, 3.20e-04]    []  \r\n",
      "Epoch 93000: train loss improved from 3.09e-03 to 3.05e-03, saving model to Model/model-93000.pdparams ...\r\n",
      "\r\n",
      "94000     [1.92e-04, 2.89e-04, 9.92e-05, 1.39e-04, 1.03e-03, 4.58e-04, 4.87e-04, 3.18e-04]    [2.62e-04, 6.85e-04, 3.26e-03, 1.34e-03, 1.03e-03, 4.58e-04, 4.87e-04, 3.18e-04]    []  \r\n",
      "Epoch 94000: train loss improved from 3.05e-03 to 3.02e-03, saving model to Model/model-94000.pdparams ...\r\n",
      "\r\n",
      "95000     [1.88e-04, 2.82e-04, 9.82e-05, 1.36e-04, 1.03e-03, 4.54e-04, 4.82e-04, 3.14e-04]    [2.63e-04, 6.38e-04, 3.31e-03, 1.27e-03, 1.03e-03, 4.54e-04, 4.82e-04, 3.14e-04]    []  \r\n",
      "Epoch 95000: train loss improved from 3.02e-03 to 2.98e-03, saving model to Model/model-95000.pdparams ...\r\n",
      "\r\n",
      "96000     [1.78e-04, 2.74e-04, 9.71e-05, 1.37e-04, 1.03e-03, 4.49e-04, 4.78e-04, 3.13e-04]    [2.66e-04, 6.30e-04, 3.33e-03, 1.30e-03, 1.03e-03, 4.49e-04, 4.78e-04, 3.13e-04]    []  \r\n",
      "Epoch 96000: train loss improved from 2.98e-03 to 2.95e-03, saving model to Model/model-96000.pdparams ...\r\n",
      "\r\n",
      "97000     [1.69e-04, 2.62e-04, 9.97e-05, 1.35e-04, 1.02e-03, 4.46e-04, 4.79e-04, 3.11e-04]    [2.84e-04, 5.98e-04, 3.38e-03, 1.44e-03, 1.02e-03, 4.46e-04, 4.79e-04, 3.11e-04]    []  \r\n",
      "Epoch 97000: train loss improved from 2.95e-03 to 2.92e-03, saving model to Model/model-97000.pdparams ...\r\n",
      "\r\n",
      "98000     [1.57e-04, 2.58e-04, 9.80e-05, 1.30e-04, 1.02e-03, 4.43e-04, 4.74e-04, 3.08e-04]    [2.61e-04, 5.85e-04, 3.48e-03, 1.36e-03, 1.02e-03, 4.43e-04, 4.74e-04, 3.08e-04]    []  \r\n",
      "Epoch 98000: train loss improved from 2.92e-03 to 2.88e-03, saving model to Model/model-98000.pdparams ...\r\n",
      "\r\n",
      "99000     [1.40e-04, 2.55e-04, 9.95e-05, 1.28e-04, 1.01e-03, 4.39e-04, 4.73e-04, 3.07e-04]    [2.54e-04, 5.55e-04, 3.87e-03, 1.29e-03, 1.01e-03, 4.39e-04, 4.73e-04, 3.07e-04]    []  \r\n",
      "Epoch 99000: train loss improved from 2.88e-03 to 2.85e-03, saving model to Model/model-99000.pdparams ...\r\n",
      "\r\n",
      "100000    [1.27e-04, 2.44e-04, 9.77e-05, 1.27e-04, 1.01e-03, 4.40e-04, 4.71e-04, 3.06e-04]    [2.88e-04, 5.97e-04, 4.36e-03, 1.53e-03, 1.01e-03, 4.40e-04, 4.71e-04, 3.06e-04]    []  \r\n",
      "Epoch 100000: train loss improved from 2.85e-03 to 2.82e-03, saving model to Model/model-100000.pdparams ...\r\n",
      "\r\n",
      "101000    [1.20e-04, 2.39e-04, 9.59e-05, 1.23e-04, 1.00e-03, 4.37e-04, 4.70e-04, 3.04e-04]    [2.97e-04, 6.09e-04, 4.71e-03, 1.65e-03, 1.00e-03, 4.37e-04, 4.70e-04, 3.04e-04]    []  \r\n",
      "Epoch 101000: train loss improved from 2.82e-03 to 2.79e-03, saving model to Model/model-101000.pdparams ...\r\n",
      "\r\n",
      "102000    [1.15e-04, 2.33e-04, 9.49e-05, 1.20e-04, 9.94e-04, 4.35e-04, 4.68e-04, 3.03e-04]    [3.12e-04, 6.43e-04, 4.70e-03, 1.86e-03, 9.94e-04, 4.35e-04, 4.68e-04, 3.03e-04]    []  \r\n",
      "Epoch 102000: train loss improved from 2.79e-03 to 2.76e-03, saving model to Model/model-102000.pdparams ...\r\n",
      "\r\n",
      "103000    [1.10e-04, 2.17e-04, 9.26e-05, 1.18e-04, 9.91e-04, 4.36e-04, 4.69e-04, 3.05e-04]    [3.48e-04, 6.64e-04, 4.35e-03, 2.12e-03, 9.91e-04, 4.36e-04, 4.69e-04, 3.05e-04]    []  \r\n",
      "Epoch 103000: train loss improved from 2.76e-03 to 2.74e-03, saving model to Model/model-103000.pdparams ...\r\n",
      "\r\n",
      "104000    [1.07e-04, 2.17e-04, 9.34e-05, 1.14e-04, 9.89e-04, 4.30e-04, 4.67e-04, 3.04e-04]    [3.72e-04, 6.65e-04, 4.44e-03, 2.22e-03, 9.89e-04, 4.30e-04, 4.67e-04, 3.04e-04]    []  \r\n",
      "Epoch 104000: train loss improved from 2.74e-03 to 2.72e-03, saving model to Model/model-104000.pdparams ...\r\n",
      "\r\n",
      "105000    [1.08e-04, 2.19e-04, 8.96e-05, 1.14e-04, 9.85e-04, 4.24e-04, 4.62e-04, 2.99e-04]    [3.09e-04, 5.94e-04, 4.33e-03, 1.90e-03, 9.85e-04, 4.24e-04, 4.62e-04, 2.99e-04]    []  \r\n",
      "Epoch 105000: train loss improved from 2.72e-03 to 2.70e-03, saving model to Model/model-105000.pdparams ...\r\n",
      "\r\n",
      "106000    [1.06e-04, 2.11e-04, 8.66e-05, 1.10e-04, 9.82e-04, 4.25e-04, 4.63e-04, 2.98e-04]    [3.12e-04, 5.61e-04, 4.45e-03, 1.89e-03, 9.82e-04, 4.25e-04, 4.63e-04, 2.98e-04]    []  \r\n",
      "Epoch 106000: train loss improved from 2.70e-03 to 2.68e-03, saving model to Model/model-106000.pdparams ...\r\n",
      "\r\n",
      "107000    [1.03e-04, 2.11e-04, 8.46e-05, 1.08e-04, 9.77e-04, 4.24e-04, 4.61e-04, 2.98e-04]    [3.27e-04, 5.29e-04, 4.65e-03, 1.96e-03, 9.77e-04, 4.24e-04, 4.61e-04, 2.98e-04]    []  \r\n",
      "Epoch 107000: train loss improved from 2.68e-03 to 2.67e-03, saving model to Model/model-107000.pdparams ...\r\n",
      "\r\n",
      "108000    [1.02e-04, 2.08e-04, 8.09e-05, 1.03e-04, 9.76e-04, 4.23e-04, 4.58e-04, 2.96e-04]    [3.31e-04, 4.77e-04, 4.97e-03, 1.89e-03, 9.76e-04, 4.23e-04, 4.58e-04, 2.96e-04]    []  \r\n",
      "Epoch 108000: train loss improved from 2.67e-03 to 2.65e-03, saving model to Model/model-108000.pdparams ...\r\n",
      "\r\n",
      "109000    [9.95e-05, 2.04e-04, 7.88e-05, 1.01e-04, 9.70e-04, 4.25e-04, 4.56e-04, 2.94e-04]    [3.41e-04, 4.55e-04, 5.20e-03, 1.92e-03, 9.70e-04, 4.25e-04, 4.56e-04, 2.94e-04]    []  \r\n",
      "Epoch 109000: train loss improved from 2.65e-03 to 2.63e-03, saving model to Model/model-109000.pdparams ...\r\n",
      "\r\n",
      "110000    [9.68e-05, 2.04e-04, 7.67e-05, 9.99e-05, 9.64e-04, 4.24e-04, 4.53e-04, 2.93e-04]    [3.86e-04, 4.42e-04, 5.50e-03, 2.16e-03, 9.64e-04, 4.24e-04, 4.53e-04, 2.93e-04]    []  \r\n",
      "Epoch 110000: train loss improved from 2.63e-03 to 2.61e-03, saving model to Model/model-110000.pdparams ...\r\n",
      "\r\n",
      "\r\n",
      "Best model at step 110000:\r\n",
      "  train loss: 2.61e-03\r\n",
      "  test loss: 1.06e-02\r\n",
      "  test metric: []\r\n",
      "\r\n",
      "'train' took 2151.887946 s\r\n",
      "\r\n",
      "Saving loss history to ./Train_history/loss.dat ...\r\n",
      "Saving training data to ./Train_history/train.dat ...\r\n",
      "Saving test data to ./Train_history/test.dat ...\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(640x480)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "Figure(1000x400)\r\n",
      "/home/aistudio\r\n"
     ]
    }
   ],
   "source": [
    "#翼型Naca0010-65，Level1网格训练\n",
    "%cd /home/aistudio/work/NACA0010-65/level1/2_PINN_training/\n",
    "!python level1_training_naca001065.py\n",
    "%cd /home/aistudio/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "execution": {
     "iopub.execute_input": "2023-11-09T03:29:16.616041Z",
     "iopub.status.busy": "2023-11-09T03:29:16.615490Z",
     "iopub.status.idle": "2023-11-09T03:29:23.080707Z",
     "shell.execute_reply": "2023-11-09T03:29:23.079485Z",
     "shell.execute_reply.started": "2023-11-09T03:29:16.616010Z"
    },
    "jupyter": {
     "outputs_hidden": true
    },
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/aistudio/work/NACA6334-18/level1/2_PINN_training\r\n",
      "Using backend: paddle\r\n",
      "Other available backends: tensorflow.compat.v1, tensorflow, pytorch, jax.\r\n",
      "paddle supports more examples now and is recommended.\r\n",
      " \r\n",
      "Set the default float type to float32\r\n",
      "Figure(1920x1440)\r\n",
      "Warning: CSGDifference.uniform_points not implemented. Use random_points instead.\r\n",
      "[4522, 4522, 4522, 4522]\r\n",
      "4522\r\n",
      "Compiling model...\r\n",
      "'compile' took 0.000235 s\r\n",
      "\r\n",
      "Warning: epochs is deprecated and will be removed in a future version. Use iterations instead.\r\n",
      "Training model...\r\n",
      "\r\n",
      "Step      Train loss                                                                          Test loss                                                                           Test metric\r\n",
      "0         [6.22e-04, 9.12e-02, 5.39e-02, 1.14e-02, 1.01e+01, 1.03e+00, 8.44e+00, 8.31e+00]    [7.19e-04, 2.70e-02, 1.59e-02, 1.51e-02, 1.01e+01, 1.03e+00, 8.44e+00, 8.31e+00]    []  \r\n",
      "^C\r\n",
      "/home/aistudio\r\n"
     ]
    }
   ],
   "source": [
    "#翼型Naca6334-18，Level1网格训练\n",
    "%cd /home/aistudio/work/NACA6334-18/level1/2_PINN_training/\n",
    "!python level1_training_naca633418.py\n",
    "%cd /home/aistudio/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 直接加载训练好的模型，并绘制预测流场"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "from work.restore_model_tools import restore_model, plot_naca001065_flow_field\n",
    "level1_naca001065_model = restore_model(model_path='/home/aistudio/work/NACA0010-65/level1/2_PINN_training/Model/model-110000.pdparams')\n",
    "plot_naca001065_flow_field(level1_naca001065_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以**Naca0010-65原网格Level1**为例，最终训练完成的模型预测的流场**速度**、**压力**和**密度**云图如下：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e951b1713e2f488a9c2d02888cb9fa91751c0742e0c049f49ba77e65177f58ce)\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/24f48f6857d043afa1a57b794b5e31de0d3e57201cf84f5fbbc73bfe4f0f64bf)\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/ea8a067a4f9f44fd8ac25c1171510a03864919ca814f4aa2a5b95129d59f1fc7)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "from work.restore_model_tools import restore_model, plot_naca633418_flow_field\n",
    "level1_naca633418_model = restore_model(model_path='/home/aistudio/work/NACA6334-18/level1/2_PINN_training/Model/model-75000.pdparams')\n",
    "plot_naca633418_flow_field(level1_naca633418_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以**Naca6334-18原网格Level1**为例，最终训练完成的模型预测的流场**速度**、**压力**和**密度**云图如下：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/f3c02baefe0d4584b1e2b80393d9512de677c352a80345a3a28f7fe2f262292c)\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/5fd2ef0d9fd943afb4247e4d595a9ef3f6bef7ab741d4377ba2ef4e54733c8e4)\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/702c5a0d6c92466cbd3336bb737f95e6acd78174b316446d9d90cdca83420160)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. 基于Edges中点的自适应网格细化算法\n",
    "\n",
    "\n",
    "网格逐级细化过程大体思路是，\n",
    "\n",
    "- 原网格先使用Fluent进行模拟，PINNs训练模型。\n",
    "- 在流场区域内布满预测点来计算PDEs的残差值大小来得到一组残差较大的网格点，在该网格点所在区域内需要进行网格的细化。\n",
    "- 之后将这组残差点加入到原有的网格点集中，生成加密后的网格。\n",
    "- 再使用Fluent进行模拟，PINNs训练模型，得到新的一组残差较大的点，加入前一个level的网格中进行加密。\n",
    "\n",
    "\n",
    "\n",
    "但在得到残差点后进行网格细化的过程中遇到了棘手的问题。若直接在需细化区域内随机采样细化的坐标点，并加入到原有的网格节点集中，再通过Delaunay三角化去生成网格的话，最终得到的网格质量将会特别差（出现极大角度钝角或者极小角度锐角的情况）。因此我们采用原网格边（Edge）的中点作为预测PINNs训练好的模型PDEs残差的点，然后再进行Delaunay三角化去生成加密后的网格，可以得到很好的网格质量和细化效果。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 基于随机采样点进行网格细化\n",
    "\n",
    "在PINNs训练完成后，在区域内**随机采样100000个点**来预测每个位置的PDEs残差，随后我们使用这些区域坐标点去计算PDEs的残差，并从中筛选出前2000个左右**对应PDEs残差最大的点**并绘制成散点图："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "from work.mesh_refinement_tools import naca001065_geom, DL_Euler_Equation_2D, get_top_k_indices, plot_naca001065_add_points\n",
    "from work.restore_model_tools import restore_model\n",
    "level1_naca001065_model = restore_model(model_path='/home/aistudio/work/NACA0010-65/level1/2_PINN_training/Model/model-110000.pdparams')\n",
    "\n",
    "X = naca001065_geom().random_points(100000)\n",
    "[f1, f2, f3, f4] = level1_naca001065_model.predict(X, operator=DL_Euler_Equation_2D)\n",
    "err_eq = np.absolute(f1) + np.absolute(f2) + np.absolute(f3) + np.absolute(f4)\n",
    "x_id = get_top_k_indices(err_eq, k=2000)\n",
    "add_points = X[x_id]\n",
    "plot_naca001065_add_points(add_points)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca6334-18\n",
    "from work.mesh_refinement_tools import naca633418_geom, DL_Euler_Equation_2D, get_top_k_indices, plot_naca633418_add_points\n",
    "X = naca633418_geom().random_points(100000)\n",
    "[f1, f2, f3, f4] = level1_naca633418_model.predict(X, operator=DL_Euler_Equation_2D)\n",
    "err_eq = np.absolute(f1) + np.absolute(f2) + np.absolute(f3) + np.absolute(f4)\n",
    "x_id = get_top_k_indices(err_eq, k=2000)\n",
    "add_points = X[x_id]\n",
    "plot_naca633418_add_points(add_points)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以注意到，残差较大的点主要集中在翼型的前缘、后缘以及激波所在的位置处。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/487485cd4e8a4e30930957e6dc14948883b7587660f948e99a5a3463e107fc81)\n",
    "\n",
    "接下来将这2000个点加入到原网格点集中，再进行约束Delauney三角化后但得到加密质量很差的网格：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/f186ad9f60a04c2f80732d3f8ead57e2c1df0f87107c45d0a7ac08685edf6a05)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 基于Edges中点的自适应网格细化\n",
    "\n",
    "为此我们提出并采用**基于三角形网格边（Edge）中点的一变四（1to4）细分网格**的算法。主要思路为：\n",
    "\n",
    "- 对每一个triangular mesh的每条边（Edge）取中点形成新的顶点，然后分别连接相邻三个顶点形成新的细分三角形。\n",
    "- 细分后的四个三角形的法线和原来三角形一致。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/18dafbc3889041afafb1d8136c84f61ebadda5819f804441a8d2df50b76c506a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于现有的翼型网格，主要流程就是计算出原网格中每条边的中点。如下图所示，黑色散点代表原网格的4511个节点，红色散点则代表原网格中每条Edge的中点（共13168个点）。\n",
    "\n",
    "网格加密即通过PINNs选出处于残差较大区域的红色中点，再添加到原网格的节点中去，重新Delauney三角化即可生成局部细化的网格。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e2bf4d6f14a243e0b21044feb6505f641901d5d0b1a44fcc96652ee0d7e3ba53)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/8524e8f8f36d4792bb1c16cbf3d1af8e848f9e708d964381a04479aa0e8e30a0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**接下来的网格细化过程展示以比赛方提供的翼型NACA0010-65的Level1原始网格为例对代码和流程进行详细说明。**\n",
    "\n",
    "两个翼型网格细化的整体代码均在各level/3_PINN_refine_mesh文件夹的PINN_guide_mesh_refinement.ipynb中"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**需要注意的是，由于带边界约束的Delauney三角化需要指定翼型边界点的索引，但此时加入的Edges中点其中也包括少量处于翼型表面位置的点。因此我们需要筛选出此时翼型边界的所有点，更新边界约束的索引。**\n",
    "因此我们在计算所有Edges的中点前，先计算翼型的中点，并加入到翼型表面原本点中，顺序排列，相当于对翼型表面进行了细化。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:30.616871Z",
     "iopub.status.busy": "2023-11-17T09:31:30.615869Z",
     "iopub.status.idle": "2023-11-17T09:31:30.625395Z",
     "shell.execute_reply": "2023-11-17T09:31:30.624485Z",
     "shell.execute_reply.started": "2023-11-17T09:31:30.616833Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from work.mesh_refinement_tools import insert_points_between_points\n",
    "#翼型Naca0010-65边界细化\n",
    "# airfoil_points1 是原翼型坐标点，形状为 (205, 2)\n",
    "# refined_airfoil_points1 是细化后的翼型坐标点，形状为 (410, 2)\n",
    "refined_airfoil_points = insert_points_between_points(airfoil_points1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后得到加密前网格的所有边的信息，使用-e选项输出Edges的信息："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:32.585667Z",
     "iopub.status.busy": "2023-11-17T09:31:32.585065Z",
     "iopub.status.idle": "2023-11-17T09:31:32.595132Z",
     "shell.execute_reply": "2023-11-17T09:31:32.594250Z",
     "shell.execute_reply.started": "2023-11-17T09:31:32.585630Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# 获取翼型Naca0010-65的level1网格的节点和边的信息\n",
    "A = dict(vertices = polygon1, segments = inner_outer_index1, holes=[[0.5, 0]])\n",
    "B = tr.triangulate(A,'pe')\n",
    "edges = B['edges']\n",
    "vertices = B['vertices']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "筛选出原网格除了翼型的边之外的所有边："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:34.364034Z",
     "iopub.status.busy": "2023-11-17T09:31:34.362857Z",
     "iopub.status.idle": "2023-11-17T09:31:34.478662Z",
     "shell.execute_reply": "2023-11-17T09:31:34.477556Z",
     "shell.execute_reply.started": "2023-11-17T09:31:34.363997Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "connected_edges = set()\n",
    "for edge_index in edges:\n",
    "    vertex1_index, vertex2_index = edge_index\n",
    "    if vertex1_index in inner_airfoil_index1 and vertex2_index in inner_airfoil_index1:\n",
    "        connected_edges.add(tuple(edge_index))  # Convert the numpy.ndarray to a tuple\n",
    "filtered_edges = [edge_index for edge_index in edges if tuple(edge_index) not in connected_edges]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "得到除翼型外所有Edges的中点："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:36.466981Z",
     "iopub.status.busy": "2023-11-17T09:31:36.465995Z",
     "iopub.status.idle": "2023-11-17T09:31:36.529905Z",
     "shell.execute_reply": "2023-11-17T09:31:36.528884Z",
     "shell.execute_reply.started": "2023-11-17T09:31:36.466944Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "inner_midpoints = []\n",
    "for edge in filtered_edges:\n",
    "    vertex1 = vertices[edge[0]]\n",
    "    vertex2 = vertices[edge[1]]\n",
    "    inner_midpoint = (vertex1 + vertex2) / 2\n",
    "    inner_midpoints.append(inner_midpoint)\n",
    "inner_midpoints = np.array(inner_midpoints)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "合并所有边的中点，翼型点在前，其他点在后："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:38.759310Z",
     "iopub.status.busy": "2023-11-17T09:31:38.758374Z",
     "iopub.status.idle": "2023-11-17T09:31:38.763377Z",
     "shell.execute_reply": "2023-11-17T09:31:38.762576Z",
     "shell.execute_reply.started": "2023-11-17T09:31:38.759275Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "all_edges_midpoints = np.vstack((refined_airfoil_points[1::2], inner_midpoints))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后我们使用原网格的Edges中点的13168个区域坐标点去计算PDEs的残差，并从中筛选出前2000个左右（约前15%）对应PDEs残差最大的点并绘制成散点图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "from work.mesh_refinement_tools import DL_Euler_Equation_2D, get_top_k_indices, plot_naca001065_add_points\n",
    "from work.restore_model_tools import restore_model\n",
    "level1_naca001065_model = restore_model(model_path='/home/aistudio/work/NACA0010-65/level1/2_PINN_training/Model/model-110000.pdparams')\n",
    "X = all_edges_midpoints\n",
    "[f1, f2, f3, f4] = level1_naca001065_model.predict(X, operator=DL_Euler_Equation_2D)\n",
    "err_eq = np.absolute(f1) + np.absolute(f2) + np.absolute(f3) + np.absolute(f4)\n",
    "x_id = get_top_k_indices(err_eq, k=2000)\n",
    "add_points = X[x_id]\n",
    "plot_naca001065_add_points(X[x_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到PDEs残差较大的Edges的中点\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/d1f59b6d97da4b67a132baa195a87bbc5b81576ece57406983c44b611def9b31)\n",
    "\n",
    "将待加入进行细化的点进行分类，分清哪些是翼型上的，哪些不是："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "更新新网格的所有点"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-11-17T09:31:46.977831Z",
     "iopub.status.busy": "2023-11-17T09:31:46.977031Z",
     "iopub.status.idle": "2023-11-17T09:31:46.989040Z",
     "shell.execute_reply": "2023-11-17T09:31:46.988064Z",
     "shell.execute_reply.started": "2023-11-17T09:31:46.977794Z"
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from work.mesh_refinement_tools import insert_midpoints\n",
    "# 处于翼型中点的那些点\n",
    "airfoil_midpoints = X[x_id[x_id < refined_airfoil_points[1::2].shape[0]]]\n",
    "# 翼型上待加密的点插到原翼型所有点里\n",
    "airfoil_newpoints = insert_midpoints(airfoil_midpoints, airfoil_points1)\n",
    "# 不包括原网格的节点的中点\n",
    "internal_midpoints =  X[x_id[x_id >= refined_airfoil_points[1::2].shape[0]]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "再以这些新的网格节点重新进行约束Delauney三角化后可生成加密后的网格，质量可得到大幅提高！"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "'''----------------------------- Delaunay triangulation for refinement mesh--------------------------'''\n",
    "polygon = np.concatenate((airfoil_newpoints, farfield_points1, internal_points1, internal_midpoints))\n",
    "# 设置约束边界的索引，翼型边界和矩形边界\n",
    "inner_airfoil_index = np.hstack((np.arange(0, airfoil_newpoints.shape[0], 1).reshape(airfoil_newpoints.shape[0], 1), \n",
    "                                  np.arange(1, airfoil_newpoints.shape[0]+1, 1).reshape(airfoil_newpoints.shape[0], 1)))\n",
    "inner_airfoil_index[airfoil_newpoints.shape[0]-1, 1]=0\n",
    "outer_farfield_index = np.hstack((np.arange(0, farfield_points1.shape[0], 1).reshape(farfield_points1.shape[0], 1), \n",
    "                                  np.arange(1, farfield_points1.shape[0] +1, 1).reshape(farfield_points1.shape[0], 1)))\n",
    "outer_farfield_index[farfield_points1.shape[0]-1, 1]=0\n",
    "outer_farfield_index = outer_farfield_index + inner_airfoil_index.shape[0]\n",
    "inner_outer_index = np.vstack((inner_airfoil_index, outer_farfield_index))\n",
    "# 三角化\n",
    "C = dict(vertices = polygon, segments = inner_outer_index, holes=[[0.5, 0]])\n",
    "D = tr.triangulate(C,'p')\n",
    "tr.comparev(plt, C, D, figsize=(60, 45))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](https://ai-studio-static-online.cdn.bcebos.com/9e4806d598d34f2c9e13223bc7b44935aba657d93e5a42e8a4d12c71358c43c0)\n",
    "\n",
    "\n",
    "可以明显地看到，**我们的自适应策略主要对翼型的前缘、后缘以及激波所在的区域进行了网格加密**，这较为符合物理直觉：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/1cf51a03dfaf41229a53c8ea27bfbbe9f1824f628962471d9419ecabc344e6de)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**至此，我们完成了一轮完整的PINNs指导的网格加密。**随后可以使用meshio将该网格转换成.vtk或者.msh文件："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "#翼型Naca0010-65\n",
    "'''----------------------------- convert midpoint_refinement mesh to vtk and msh --------------------------'''\n",
    "import meshio\n",
    "\n",
    "level = 1\n",
    "\n",
    "# 获取三角剖分结果的顶点和单元格信息\n",
    "vertices = D['vertices']\n",
    "cells = [(\"triangle\", D['triangles'])]\n",
    "\n",
    "# 创建meshio的网格对象\n",
    "mesh = meshio.Mesh(points=vertices, cells=cells)\n",
    "\n",
    "# 定义保存的文件名\n",
    "vtkFileName = 'level' + str(level+1) + '_midpoints_refinement_mesh.vtk'\n",
    "# 将网格对象保存为VTK文件\n",
    "meshio.write(vtkFileName, mesh, file_format=\"vtk\", binary=False)\n",
    "\n",
    "# 使用meshio读取VTK文件\n",
    "mesh = meshio.read(vtkFileName)\n",
    "\n",
    "# 将网格对象保存为Gmsh支持的格式（例如.msh文件）\n",
    "mshFileName = 'level' + str(level+1) + '_midpoints_refinement_mesh.msh'\n",
    "meshio.write(mshFileName, mesh, file_format='gmsh22',binary=False)\n",
    "print(\"转换后的文件已保存为:\", mshFileName)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之后的不同网格层级之间的加密均采用此细化策略。所有代码和文件均在work文件夹中，此处不详细展开。**注意，代码生成的网格.msh文件仍然需要使用gmsh将其转化为.bdf文件，然后再导入Workbench中设置边界类型供Fluent计算**。\n",
    "\n",
    "\n",
    "\n",
    "综上，按照以上思路和流程，我们总共对翼型NACA0010-65的原网格Level1 → Level2 → Level3 → Level4 → Level5 **进行了四轮网格加密**，最终通过升阻力的对比与分析，认为Level4的网格为最优网格。\n",
    "\n",
    "我们总共对翼型NACA6334-18的原网格Level1 → Level2 → Level3 → Level4 **进行了三轮网格加密**，最终通过升阻力的对比与分析，认为Level2的网格为最优网格。\n",
    "\n",
    "**详细的升阻力和翼型表面的压力分布将在下一节中进行展示**。\n",
    "\n",
    "\n",
    "\n",
    "#### 1. 以下是NACA0010-65不同的优化Level所对应的翼型附近的网格：\n",
    "\n",
    "#### Level1（4,511 Points，8,657 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/5c8ce9e9f29749adb0e10b9bde6141941a8983815ffb4047a5f6e61188274dc4)\n",
    "\n",
    "#### Level2（6,397 Points，12,308 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/2532c1fb250841da92d89bb7a8a2995742f89149ca7441919ffb0faf674c5f23)\n",
    "\n",
    "#### Level3（8,293 Points，16,042 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/91a3817545e244fab5df11de3acd3ece4beea5956b4d44f78c4adecbdedf29bc)\n",
    "\n",
    "#### Level4（10,403 Points，20,208 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/00c62a023d5c4fc2ba5d37abb589270cdc06163fd1374194b76fa433acfe740d)\n",
    "\n",
    "#### Level5（13,421 Points，26,243 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/466a71c415c64f3aa0b506a6ef24f86ce23ec03cd9994df9baedc02291077896)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. 以下是NACA6334-18的优化Level所对应的翼型附近的网格：\n",
    "\n",
    "#### Level1（4,511 Points，8,657 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/7bc9a48a500b4f27b656beac489b5d6945ce95bbe7664a958c0925b495279aa1)\n",
    "\n",
    "#### Level2（6,397 Points，12,308 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e5e2ba284d0c40458dee2eaa8e922921d3e3dc7574d3477d91288d7b9f79fd2c)\n",
    "\n",
    "#### Level3（8,293 Points，16,042 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/7c56173b1abb46519c8d14f889388e235d24c4d983e848bf81491121e362aeb0)\n",
    "\n",
    "#### Level4（10,403 Points，20,208 Cells）\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/05b70d82642f43ce9fb1c0bf34ff92977c6396d4c6d74d6fa438448d631a6dfd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 三. 验证最优网格\n",
    "### 1.Naca0010-65\n",
    "#### Fluent自适应网格加密设置\n",
    "\n",
    "考虑到结果的准确性，所以**采用了Fluent自带的自适应网格加密方法作为参考和对比**。自适应网格加密方法可以手动设置，因为这个案例涉及可压缩气体，对密度的变化比较敏感，可以设置基于密度梯度的变化来自适应网格加密，选择预定义标准——Aerodynamics——Shock Indicator——Density-based，对细化标准以及粗化标准也可以自己定义，需要注意设置好合适的频率（迭代），结果收敛后再进行下一次的自适应网格加密。\n",
    "\n",
    "\n",
    "对原网格一共进行了**5次**自适应网格加密，最后得到了664805个网格，自适应加密后的网格可以看到对激波的位置加密情况比较明显，翼型周围的网格也进行了不同程度的加密，并且也有一定地方进行了粗化。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/3b029a2466f444d29af0685356490f859900867348fe47039e5bec43ecfa898a)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/fe44fc461c0b40c3bf97fc4da29720eb495bf264d6794523beb561eecd5d269b)\n",
    "\n",
    "\n",
    "对每次自适应加密后的网格按照与原网格一样的边界条件进行模拟，依次得到了每次加密后的升阻力结果，将其与基于Edges中点的网格细化算法加密后的升阻力结果进行对比。如下图所示，结果表明我们的基于Edges中点的网格细化算法以及fluent自适应加密算法，均会使得升阻力的值随着网格加密即网格数量的增加趋向于收敛。值得一提的是，通过与Fluent自适应算法的对比，**我们的网格细化算法不但最终可以达到与Fluent相同的精度，也可以大幅提高效率，仅用少数网格即可达到Fluent自适应加密后的大量网格的效果**。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/bbc03e6109ff4e3eb8a084dcbe1cdf44d4f21185b2314b8499c7927b833070d1)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/9e9813b9820646c6a1791ece9005b7cb18e80ae1297f4ce6a257effa5dc55979)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.1 网格结果定量分析\n",
    "\n",
    "接下来我们对最优网格进行敏感性分析，即将最优网格（Level4）进行网格数目的加倍和减半，再对比CFD结果的变化情况来判定网格的优劣。根据与Fluent自适应网格加密的结果对比，得到了基于Edges中点的网格细化算法下的最优网格，并对最优网格进行网格数增加与删减处理，最终得到三种网格，分别定义为网格N0（减少），N1（最优网格）和N2（增加）。三种网格分别进行了局部放大，可以看到对激波以及翼型周围加密比较明显，而且网格加密的质量也较好。\n",
    "\n",
    "##### 网格N0：\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/d803bf47cece4102ad9b4d7873b15345996323e5abe94cf39ed981fb792a288f)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/392eac6a7e7e4a96b164c31c717a6b9c9699d8ac0a784d39afde426eafd71802)\n",
    "\n",
    "#### 网格N1：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/02c909ba7f7e4c7e9357898acefc84f8c26bd4141a044f07b66a5713efc5b3c2)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/e1f641ad3fb24715b773c4931cdc4a3e3af9c80a83274e6b80a50c1a0d9dda99)\n",
    "\n",
    "#### 网格N2：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/7478774948d848dfb4bce411d57e8ca751f472e7fbce49d7b24078f1d1408e5f)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/08c1442047964987a190ba87fe189cc33c9c419c8fb44582ab4e11b5b016bcf5)\n",
    "\n",
    "我们分别计算了这三种网格的升阻力情况，得到的结果显示网格N1与网格数增加接近三倍的网格N2升阻力大小几乎相同，差异<1%；而N1与网格数减小两倍多的网格N0的升阻力大小对比差异在10%左右，证明N1为最优网格。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/f6027733d92846ae96d0d8f35e036969223ed0eac5614710876a3415d1aae041)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/170c3e45dd1e4db28ef4ce14c5bf2f287dae41186c0d4dd69c89d32972c57654)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1.2 网格结果定性分析\n",
    "\n",
    "同时可以得到翼型表面的压力分布，将得到的网格N1（最优网格）与网格N0（减少）和网格N2（增加）的压力分布进行定性对比，会发现**N1与N2的压力分布区别很小**，只有少数点存在差别；**而N1与N0的压力分布区别整体都比较明显，尤其是激波的位置会更加明显**。若将三种网格放在一起对比会更明显看到网格N0与网格N1, N2的压力分布的区别。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/60c032b104d04f23ab6bef1d42ede7f0f988cf32312e4ea19a7b48f020ce9fbb)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/442898f8750b4be1a01e9876bf60635c0fb7cbf8cacd419a962be712bcbb7413)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/8959bb2e6b4942ccaeab5db310dfe0b7eb4db0b08ee74561bf6b3cf3b968ed33)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.Naca6334-18\n",
    "#### Fluent自适应网格加密设置\n",
    "\n",
    "对原网格一共进行了**5次**自适应网格加密，最后得到了157095个网格。和Naca0010-65翼型操作一致，对每次自适应加密后的网格按照与原网格一样的边界条件进行模拟，依次得到了每次加密后的升阻力结果，将其与基于Edges中点的网格细化算法加密后的升阻力结果进行对比。如下图所示，结果表明我们的基于Edges中点的网格细化算法以及fluent自适应加密算法，均会使得升阻力的值随着网格加密即网格数量的增加趋向于收敛。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/dd44a48e4c0d43eebaa7a356ccf64d9093b6af7c191246fb99ea7da2b45b01d3)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/a90ea63d611c4d9f808894843548be368ffd2a85484f408eb3b7eb4329d600a5)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.1 网格结果定量分析\n",
    "\n",
    "接下来我们对最优网格进行敏感性分析，即将最优网格（Level2）进行网格数目的增加和删减，再对比CFD结果的变化情况来判定网格的优劣。根据与Fluent自适应网格加密的结果对比，得到了基于Edges中点的网格细化算法下的最优网格，并对最优网格进行网格数增加与删减处理，最终得到三种网格，分别定义为网格N0（减少），N1（最优网格）和N2（增加）。三种网格分别进行了局部放大，可以看到对激波以及翼型周围加密比较明显，而且网格加密的质量也较好。\n",
    "\n",
    "#### 网格N0：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/80be97faaf1d4d99bda874c627772bdf0c7bf312c6014e20ae17f79f89769381)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/7bc9a48a500b4f27b656beac489b5d6945ce95bbe7664a958c0925b495279aa1)\n",
    "\n",
    "#### 网格N1：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/ff28c2ff2db04e9886a7570aeed52eb876dcb2e1ac694714bd6bf31cab751ede)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/07272b8efe9a4db0ac07f6b38282021693f1ea768b454df78c24d604580abce9)\n",
    "\n",
    "\n",
    "#### 网格N2：\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/f2e75b04ae2542479a50b95973b787b1c0c11cdc55c84d8f91fd68e4cd94eb98)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/1e2a6962f86e464f925c06e27f5618d71a7ae1a3b8524994b49267e5d613fab1)\n",
    "\n",
    "\n",
    "我们分别计算了这三种网格的升阻力情况，得到的结果显示网格N1与网格数增加接近两倍的网格N2升阻力大小几乎相同，差异<0.1%；而N1与网格数减小接近两倍的网格N0的升阻力大小对比差异在2%左右，这个应该是由于初始网格就比较好，使得差异变化相对较小些，不过也能证明N1为最优网格。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/19fb6fbe260147a6a255ea1d300b65d94ff52b0bd46242658116dd7dac374548)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/d859eb75d4f544fbb1eb9312acc47d9780fd3f890ff64faea6f779356a49c6e6)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.2 网格结果定性分析\n",
    "\n",
    "同时可以得到翼型表面的压力分布，将得到的网格N1（最优网格）与网格N0（减少）和网格N2（增加）的压力分布进行定性对比，会发现**N1与N2的压力分布区别很小**，只有少数点存在差别；**而N1与N0的压力分布区别在前缘和尾缘以及激波位置都比较明显**。若将三种网格放在一起对比会更明显看到网格N0与网格N1, N2的压力分布的区别。\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/21c608e10d864f27b52de2449288d5f7b4f8b1ee57d5413ba41950bafa026864)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/9711b8459363442c9429d818924ff92e49369a25bfcc4dfeaaea9e0f32a916f1)\n",
    "\n",
    "![](https://ai-studio-static-online.cdn.bcebos.com/fc249636a10c46489f7c5bd20b186a1cba6de6c9b03544b7b551889f8579e50b)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "请点击[此处](https://ai.baidu.com/docs#/AIStudio_Project_Notebook/a38e5576)查看本环境基本用法.  <br>\n",
    "Please click [here ](https://ai.baidu.com/docs#/AIStudio_Project_Notebook/a38e5576) for more detailed instructions. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "py35-paddle1.2.0"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
